home *** CD-ROM | disk | FTP | other *** search
/ Aminet 2 / Aminet AMIGA CDROM (1994)(Walnut Creek)[Feb 1994][W.O. 44790-1].iso / Aminet / util / misc / Fudgit233.lha / Source / src / fft.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-14  |  6.0 KB  |  234 lines
  1. #include <math.h>
  2. #include <stdio.h>
  3.  
  4. #include "symbol.h"
  5. #include "code.h"
  6. #include "math.tab.h"
  7. #include "fudgit.h"
  8. #include "setshow.h"
  9. #include "dalloca.h"
  10. #include "head.h"
  11.  
  12. static void fft(double *data, int nn, int isign);
  13. static void merge(double *real, double *ima, double *comp, int nn);
  14. static void split(double *comp, double *real, double *ima, int nn);
  15. static int checkbit(unsigned int data);
  16. static int fftsmooth(double *fromvec, double *smvec, int ndata, float sigma);
  17.  
  18. /* called with 4 arguments (argc=5) */
  19. int Ft_run_fft(int argc, char **argv, int dir, int ndata)
  20. {
  21.     double *vectors[4], *complex;
  22.     Symbol *sym;
  23.     int i;
  24.     extern Symbol *Ft_lookup(char *s), *Ft_install(char *s, int t, int size);
  25.     int Ft_varcpy(char *to, char *from);
  26.  
  27.     if (ndata == 0) {
  28.         fputs("(inv)fft: No data!\n", stderr);
  29.         return(ERRR);
  30.     }
  31.     if (checkbit(ndata) != 1) {
  32.         fputs("(inv)fft: Fast fourier transform requires 2^n data points.\n",
  33.         stderr);
  34.         return(ERRR);
  35.     }
  36.     ADVECTOR(complex, 1, 2*ndata, "(inv)fft", return);
  37.     for (i=1; i<=4; i++) {
  38.         if ((sym = Ft_lookup(argv[i])) == 0) { /* get all vectors */
  39.             if (i<=2) {
  40.                 fprintf(stderr, "(inv)fft: %s: Unknown vector.\n", argv[i]);
  41.                 return(ERRR);
  42.             }
  43.             if (Ft_varcpy(NULL, argv[i]) != VEC) {
  44.                 fprintf(stderr,
  45.                 "(inv)fft: %s: Not a legal vector name.\n", argv[i]);
  46.                 return(ERRR);
  47.             }
  48.             sym = Ft_install(argv[i], VEC, Ft_Samples);
  49.         }
  50.         vectors[i-1] = sym->u.vec;
  51.     }
  52.     /* make it complex */
  53.     merge(vectors[0], vectors[1], complex, ndata);
  54.     fft(complex, ndata, dir);  /* take the fourier transform  */
  55.     /* split in two vectors  */
  56.     split(complex, vectors[2], vectors[3], ndata);
  57.     return(0);
  58. }
  59.  
  60. /* called with 3 arguments (argc=4) */
  61. int Ft_run_smooth(int argc, char **argv, int ndata)
  62. {
  63.     double *vectors[2];
  64.     Symbol *sym;
  65.     int i;
  66.     float sigma;
  67.     extern Symbol *Ft_lookup(char *s), *Ft_install(char *s, int t, int size);
  68.     int Ft_varcpy(char *to, char *from);
  69.  
  70.     if (ndata == 0) {
  71.         fputs("smooth: No data!\n", stderr);
  72.         return(ERRR);
  73.     }
  74.     if (sscanf(argv[1], "%f", &sigma) != 1) {
  75.         fputs("Error: smooth: Cannot read smoothing factor.\n", stderr);
  76.         return(ERRR);
  77.     }
  78.     for (i=2; i<=3; i++) {
  79.         if ((sym = Ft_lookup(argv[i])) == 0) { /* get all vectors */
  80.             if (i==2) {
  81.                 fprintf(stderr,
  82.                 "Error: smooth: %s: Unknown vector.\n", argv[i]);
  83.                 return(ERRR);
  84.             }
  85.             if (Ft_varcpy(NULL, argv[i]) != VEC) {
  86.                 fprintf(stderr,
  87.                 "Error: smooth: %s: Not a legal vector name.\n", argv[i]);
  88.                 return(ERRR);
  89.             }
  90.             sym = Ft_install(argv[i], VEC, Ft_Samples);
  91.         }
  92.         vectors[i-2] = sym->u.vec;
  93.     }
  94.     return(fftsmooth(vectors[0], vectors[1], ndata, sigma));
  95. }
  96.  
  97. static int checkbit(unsigned int data)
  98. {
  99.     int bitsum = 0;
  100.     int i;
  101.  
  102.     for (i=0;i<8*sizeof(int);i++) {
  103.         bitsum += data & 01;
  104.         data >>= 1;
  105.     }
  106.     return(bitsum);
  107. }
  108.  
  109. static void split(double *comp, double *real, double *ima, int nn)
  110. {
  111.     int i, j;
  112.     double norm, sqrt(double);
  113.  
  114.     norm = sqrt((double)nn);
  115.     for (i=1,j=1;j<=nn;i+=2,j++) {
  116.            real[j] = comp[i]/norm;
  117.            ima[j] = comp[i+1]/norm;
  118.     }
  119.     return;
  120. }
  121.  
  122. static void merge(double *real, double *ima, double *comp, int nn)
  123. {
  124.     int i, j;
  125.  
  126.     if (ima) {
  127.         for (i=j=1;j<=nn;i+=2,j++) { /* merge, i always odd */
  128.             comp[i] = real[j];  /* real  */
  129.             comp[i+1] = ima[j];  /* imaginary  */
  130.         }
  131.     }
  132.     else {
  133.         for (i=j=1;j<=nn;i+=2,j++) { /* merge, i always odd */
  134.             comp[i] = real[j];  /* real  */
  135.             comp[i+1] = 0.0;  /* imaginary  */
  136.         }
  137.     }
  138.     return;
  139. }
  140.  
  141. #define TWO_PI  6.2831853071795864769252868
  142.  
  143. static void fft(double *data, int nn, int isign)
  144. {
  145.     int n, mmax, m, j, istep, i;
  146.     double wtemp, wr, wpr, wpi, wi, theta;
  147.     double rtemp, itemp;
  148.  
  149.     n = nn<<1;
  150.     j = 1;
  151.     for (i=1;i<n;i+=2) {
  152.         if (j > i) {
  153.             rtemp = data[j];
  154.             data[j] = data[i];
  155.             data[i] = rtemp;
  156.             itemp = data[j+1];
  157.             data[j+1] = data[i+1];
  158.             data[i+1] = itemp;
  159.         }
  160.         m = n>>1;
  161.         while (m >= 2 && j > m) {
  162.             j -= m;
  163.             m >>= 1;
  164.         }
  165.         j += m;
  166.     }
  167.     mmax = 2;
  168.     while (n > mmax) {
  169.         istep = mmax<<1;
  170.         theta = TWO_PI/(isign*mmax);
  171.         wtemp = sin(0.5*theta);
  172.         wpr = -2.0*wtemp*wtemp;
  173.         wpi = sin(theta);
  174.         wr = 1.0;
  175.         wi = 0.0;
  176.         for (m=1;m<mmax;m+=2) {
  177.             for (i=m;i<=n;i+=istep) {
  178.                 j = i+mmax;
  179.                 rtemp = wr*data[j] - wi*data[j+1];
  180.                 itemp = wr*data[j+1] + wi*data[j];
  181.                 data[j] = data[i]-rtemp;
  182.                 data[j+1] = data[i+1]-itemp;
  183.                 data[i] += rtemp;
  184.                 data[i+1] += itemp;
  185.             }
  186.             wr = (wtemp=wr) * wpr - wi * wpi + wr;
  187.             wi = wi * wpr + wtemp * wpi+ wi;
  188.         }
  189.         mmax = istep;
  190.     }
  191.     return;
  192. }
  193.  
  194. static int fftsmooth(double *fromvec, double *smvec, int ndata, float sigma)
  195. {
  196.     int fmax, num, k, kk, f, last, vecsize=2;
  197.     double yyn, yy1, rn1, fac;
  198.     double *comp;
  199.  
  200.     k = ndata << 1;
  201.     while (vecsize < k)
  202.         vecsize <<= 1;
  203.     ADVECTOR(comp, 1, vecsize, "smooth", return);
  204.     merge(fromvec, NULL, comp, ndata);
  205.     num = vecsize >> 1;      /* num = number of complex elements */
  206.     fmax = vecsize >> 2;    /* maximum frequency  */
  207.     sigma *= fmax;
  208.     sigma *= -sigma;
  209.     last = 2*ndata-1;
  210.     yy1 = comp[1];
  211.     yyn = comp[last];
  212.     rn1 = 1.0/(ndata-1);
  213.     for (kk=f=1;f<=ndata;f++,kk+=2) {
  214.         comp[kk] += (-rn1*(yy1*(ndata-f) + yyn*(f-1)));
  215.     }
  216.     for (f=last+2;f<=vecsize;f++) {
  217.         comp[f] = 0.0;   /* zero pad both real and ima */
  218.     }
  219.     fft(comp, num, 1);
  220.     comp[1] /= num;  /* normalize: should be real (comp[2] == 0) */
  221.     for (f=1;f<=fmax;f++) {
  222.         k = 2*f+1;                    /* positive frequencies */
  223.         kk = vecsize-k+2;                /* negative frequencies */
  224.         fac = exp(f*f/sigma)/num;     /* Gaussian window */
  225.         comp[kk] = (comp[k] *= fac);  /* real */
  226.         comp[kk+1] = -(comp[k+1] *= fac);  /* ima */
  227.     }
  228.     fft(comp, num, -1);
  229.     for (kk=f=1;f<=ndata;f++,kk+=2) {
  230.         smvec[f] = comp[kk] + (rn1*(yy1*(ndata-f) + yyn*(f-1)));
  231.     }
  232.     return(0);
  233. }
  234.